ApoE4 dose effects on serum metabolome in Alzheimer’s Disease

Global test in AD group

Performing a Global Test on the serum metabolites of AD patients, correcting for sex (Ho: E4 status has no effect on metabolite levels, Ha: it has an effect), yields a significant difference (p = 0.029 ) between E4 carriers and non-carriers. The presence of ApoE4, correcting for sex seems to have an effect on the metabolites shown in the plot below. The most significantly affected metabolites are triglycerides, diglycerides (positively) and γ-aminobutyric acid (negatively).

load("data/clinical.Rdata")
clinical$V_MMSE <- NULL
df[,1:230] <- scale(df[,1:230])
AD <- subset(df, Diagnosis == "Probable AD")
AD$sex <- as.numeric(AD$sex) -1
# Testing for effects
gt.b <- globaltest::gt(E4 ~ 1, E4 ~ .-E4dose -Diagnosis -target, data = AD)
gt.b

  p-value Statistic Expected Std.dev #Cov
1  0.0293      1.35    0.794   0.242  231

globaltest::covariates(gt.b)

Performing a Global Test on ApoE4 dose shows distinct effects depeding on the number of ApoE4 alleles (p = 0.199).

# Multinomial outcome
AD$E4dose <- as.factor(AD$E4dose)
gt.m <- globaltest::gt(E4dose ~ 1, E4dose ~ .-E4 -Diagnosis -target, data = AD)
gt.m

  p-value Statistic Expected Std.dev #Cov
1  0.0199      1.22      0.8   0.169  231

globaltest::covariates(gt.m)

Nested Linear Models

Metabolite-level models to test for ApoE4 dose effects. The F-tests yielded p-values < 0.05, however none of them survived FDR correction.

nested <- function(Y, x, ...) {
  ### Analysis of Variance (ANOVA)
  x <- as.factor(x)
  df <- cbind(Y, x, ...)
  ncol <- ncol(Y)
  covariates <- names(...)
  F_tests <- furrr::future_map(df[, 1:ncol], ~ {
    frm <- as.formula(paste0(".x ~", paste(covariates, collapse = "+")))
    rm <- lm(frm, data = df)
    ffm <- as.formula(paste0(".x ~ x +", paste(covariates, collapse = "+")))
    fm <- lm(ffm, data = df)
    anova(rm, fm)
  })
  # Correction for Multiple Testing
  # Create a list to store the p-values
  p_values <- list(1:ncol)
  # Extract the p-values of the F-tests from the anova summaries list and store them in p_values
  for (i in 1:ncol) {
    p_values[[i]] <- F_tests[[i]][["Pr(>F)"]][[2]]
  }
  # Coerce p_values to dataframe and transpose it
  p_values <- t(data.frame(p_values))

  # Set as row names the metabolites
  p_values <- data.frame(p_values, row.names = colnames(Y))

  # Calculate the FDR-adjusted p-values
  p_values$p_adj <- p.adjust(p_values[, ], method = "fdr", n=230 )
  
  # Filter out the non-significant (a=0.05) FDR-adjusted p-values
  sig <- as.data.frame(dplyr::filter(p_values, p_values < 0.05))
  sig <- sig[order(sig$p_adj), ]
  cat("Significant F-tests: \n")
  print(sig)
  cat("\nDose effects on metabolites:")
  knitr::kable(sig)
  summaries <- furrr::future_map(df[, rownames(sig)], ~ {
    f <- as.formula(paste0(".x ~ x +", paste(covariates, collapse = "+")))
    mdl <- lm(f, data = df)
    summary(mdl)$coefficients[1:length(table(x)), ]
  })
  summaries <- lapply(summaries, function(x) {
    x <- x[, -c(2, 3)]
    x <- t(x)
  })
  summaries <- plyr::ldply(summaries, id = 1:2)
  return(knitr::kable(summaries, table.attr = "style='width:30%;'"))
}

SCD and AD

The metabolites most affected by ApoE4 dose are D-L-3-aminoisobutyric.acid, histamine, L-serine, ceramide d18:1/24:0, lysophosphatidylcholine 18:1, phosphatidylcholine 34:4, 36:6 and 40:4, sphingomyelin d18:1/22:0 and d18:1/24:0, LpH.PGE2, cLpH.PGA2, Lyso-phosphatidic acid 20:5 and taurine.

Y <- df[, 1:230]
df$E4dose <- as.numeric(df$E4dose) -1
df$E4dose[df$E4dose == 0] <- "No ApoE4"
df$E4dose[df$E4dose == 1] <- "1 ApoE4"
df$E4dose[df$E4dose == 2] <- "2 ApoE4"
E4dose <- df$E4dose <- as.factor(df$E4dose)
## Testing for effects of ApoE4 dose
nested(Y, relevel(E4dose, ref= "No ApoE4"), clinical)

Significant F-tests: 
                                p_values     p_adj
Am.DL.3.aminoisobutyric.acid 0.019274304 0.6617562
Am.Histamine                 0.010564489 0.6617562
Am.L.Serine                  0.029824862 0.6617562
Lip.Cer.d18.1.24.0.          0.021905538 0.6617562
Lip.LPC.O.18.1.              0.031940407 0.6617562
Lip.PC.34.4.                 0.034526409 0.6617562
Lip.PC.40.4.                 0.030049244 0.6617562
Lip.SM.d18.1.22.0.           0.028799960 0.6617562
Lip.SM.d18.1.24.0.           0.032508286 0.6617562
OS.LpH.PGE2                  0.033150673 0.6617562
OS.cLpH.PGA2                 0.031558866 0.6617562
OS.HpH.LPA.C20.5             0.009384457 0.6617562
Am.Taurine                   0.048562691 0.7013423
Lip.PC.O.36.6.               0.045023526 0.7013423

Dose effects on metabolites:

AD group only

Among AD patients, ApoE4 dose seems to have an effect on several triglycerides, diglycerides, L-glutamic acid, L-lysine, 2- and 3-hydroxybutyric acid, glutamic acid, 3-phosphoglyceric acid, phosphatidylcholin, sphyngomyelin, HpH.PAF.C17, phosphoethanolamine, imino-diacetate.

AD <- subset(df, Diagnosis == "Probable AD")
ADY <- AD[, 1:230]
ADE4dose <- AD$E4dose
ADclinical <- clinical[df$Diagnosis == "Probable AD",]

## Testing for effects of ApoE4 dose
nested(ADY, relevel(ADE4dose, ref= "No ApoE4"), ADclinical)

Significant F-tests: 
                                    p_values     p_adj
Lip.TG.52.3.                     0.002128748 0.2496400
Lip.TG.52.4.                     0.002585849 0.2496400
Lip.TG.54.5.                     0.004966959 0.2496400
Lip.TG.54.6.                     0.005426957 0.2496400
Lip.TG.56.6.                     0.003464999 0.2496400
Lip.TG.54.4.                     0.007750811 0.2845736
Lip.DG.36.3.                     0.008660937 0.2845736
Am.L.Glutamic.acid               0.011127623 0.3141256
Am.L.Lysine                      0.028681036 0.3141256
OA.OA01...2.hydroxybutyric.acid  0.020666789 0.3141256
OA.OA04...Glutamic.Acid          0.025302325 0.3141256
OA.OA18...3.Phosphoglyceric.acid 0.024813628 0.3141256
Lip.TG.50.2.                     0.027348399 0.3141256
Lip.TG.50.3.                     0.024166483 0.3141256
Lip.TG.52.2.                     0.023188871 0.3141256
Lip.TG.52.5.                     0.021890331 0.3141256
Lip.TG.56.7.                     0.017149208 0.3141256
Lip.TG.56.8.                     0.023831512 0.3141256
Lip.Cer.d18.1.24.0.              0.028242776 0.3141256
Lip.PC.O.34.3.                   0.020831975 0.3141256
OS.HpH.LPA.C18.3                 0.015373446 0.3141256
Lip.TG.51.3.                     0.031204970 0.3262338
Lip.TG.50.1.                     0.036088644 0.3451584
Lip.SM.d18.1.22.0.               0.036844964 0.3451584
OS.HpH.PAF.C16.0                 0.037517223 0.3451584
OA.OA17...3.Hydroxybutyric.acid  0.039302817 0.3476788
Lip.TG.54.3.                     0.041634484 0.3546641
Am.O.Phosphoethanolamine         0.043904743 0.3606461
OA.OA23...Iminodiacetate         0.047339603 0.3754520

Dose effects on metabolites:

SCD group only

Among individuals with SCD, lipid metabolites were not affected as much as in the AD group, with only LPC and PC showing a difference. Aminoacids (glycine, taurine, L-serine, serotonine, L-glutamic acid, L-glutamic acid, phosphoethanolamine) and oxidative stress compounds were mostly affected in this group.

SCD <- df[df$Diagnosis == "Subjectieve klachten",]
SCDY <- SCD[, 1:230]
SCDE4dose <- SCD$E4dose
SCDclinical <- clinical[df$Diagnosis == "Subjectieve klachten",]

## Testing for effects of ApoE4 dose
nested(SCDY, relevel(SCDE4dose, ref= "No ApoE4"), SCDclinical)

Significant F-tests: 
                            p_values     p_adj
Am.Glycine               0.003514740 0.3473546
Am.Taurine               0.002948325 0.3473546
OS.cLpH.PGA2             0.005011479 0.3473546
OS.HpH.LPA.C18.1         0.006040949 0.3473546
OS.cHpH.aLPA.C16.1       0.007960965 0.3662044
Am.Ethanolamine          0.013509392 0.4438800
OA.OA11...Succinic.acid  0.013146702 0.4438800
Am.L.Serine              0.017478101 0.4580221
Am.Serotonine            0.017922604 0.4580221
Am.O.Phosphoethanolamine 0.021914277 0.5040284
OS.cLpH.PGF2a            0.024487586 0.5120132
OS.LpH.PGE2              0.028774008 0.5515018
Am.L.Glutamic.acid       0.037526747 0.6165108
Lip.LPC.18.1.            0.035209617 0.6165108
OS.HpH.S.1.P.C18.1       0.040427658 0.6198908
Am.L.Aspartic.acid       0.049024437 0.6411695
Lip.PC.36.3.             0.045474810 0.6411695

Dose effects on metabolites:

Fitting a benchmark model: clinical characteristics only

benchmark <- multifit(
  X = clin_dummy,
  y = y,
  model = "multinom",
  ctrl = ctrl,
  grid = expand.grid(decay = 0),
  trace = F
)

evaluate(benchmark)

Multi-class area under the curve: 0.8143

Confusion Matrix and Statistics

          Reference
Prediction   AD ADE4  SCD SCDE4
     AD    4756 2251    2   102
     ADE4  3736  904   31   119
     SCD     34  199 1892  4158
     SCDE4   74   46 2075  4321

Overall Statistics
                                          
               Accuracy : 0.4807          
                 95% CI : (0.4744, 0.4869)
    No Information Rate : 0.3522          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.2972          
                                          
 Mcnemar's Test P-Value : < 2.2e-16       

Statistics by Class:

                     Class: AD Class: ADE4 Class: SCD Class: SCDE4
Sensitivity             0.5530      0.2659     0.4730       0.4967
Specificity             0.8537      0.8176     0.7879       0.8628
Pos Pred Value          0.6688      0.1887     0.3011       0.6631
Neg Pred Value          0.7815      0.8746     0.8855       0.7592
Precision               0.6688      0.1887     0.3011       0.6631
Recall                  0.5530      0.2659     0.4730       0.4967
F1                      0.6054      0.2208     0.3680       0.5680
Prevalence              0.3482      0.1377     0.1619       0.3522
Detection Rate          0.1926      0.0366     0.0766       0.1749
Detection Prevalence    0.2879      0.1939     0.2544       0.2638
Balanced Accuracy       0.7034      0.5417     0.6304       0.6797

Regularized Multinomial Regression adding 230 metabolites

mlr <- multifit(
  X = Xclin,
  y = y,
  model = "multinom",
  ctrl = ctrl,
  grid = expand.grid(decay = 10),
  trace = F
)

evaluate(mlr)

Multi-class area under the curve: 0.8338

Confusion Matrix and Statistics

          Reference
Prediction   AD ADE4  SCD SCDE4
     AD    5764 2368    0     5
     ADE4  2829  966    1    42
     SCD      4    6 1659  2836
     SCDE4    3   60 2340  5817

Overall Statistics
                                          
               Accuracy : 0.5751          
                 95% CI : (0.5689, 0.5813)
    No Information Rate : 0.3522          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.4076          
                                          
 Mcnemar's Test P-Value : < 2.2e-16       

Statistics by Class:

                     Class: AD Class: ADE4 Class: SCD Class: SCDE4
Sensitivity             0.6702     0.28412    0.41475       0.6686
Specificity             0.8526     0.86516    0.86251       0.8498
Pos Pred Value          0.7084     0.25169    0.36826       0.7077
Neg Pred Value          0.8288     0.88333    0.88408       0.8251
Precision               0.7084     0.25169    0.36826       0.7077
Recall                  0.6702     0.28412    0.41475       0.6686
F1                      0.6888     0.26692    0.39012       0.6876
Prevalence              0.3482     0.13765    0.16194       0.3522
Detection Rate          0.2334     0.03911    0.06717       0.2355
Detection Prevalence    0.3294     0.15538    0.18239       0.3328
Balanced Accuracy       0.7614     0.57464    0.63863       0.7592

Projection to Latent Factors

project <- function(X, m, seed) {
  set.seed(seed)
  X <- scale(as.matrix(df[,1:230]))
  cov <- cor(X)

  # Find redundant features
  filter <- RF(cov)

  # Filter out redundant features
  filtered <- subSet(X, filter)

  # Regularized correlation matrix estimation
  M <- regcor(filtered)

  # Get the regularized correlation matrix of the filtered dataset
  R <- M$optCor

  mlfa <- mlFA(R, m = 6)
  thomson <- facScore(filtered, mlfa$Loadings, mlfa$Uniqueness)
  return(thomson)
}

thomson <- project(X, 6, seed = 1234)
X <- cbind(clin_dummy, thomson)

Multinomial Logistic Regression adding 6 factors

mlrf <- multifit(
  X = X,
  y = y,
  model = "multinom",
  ctrl = ctrl,
  trace = FALSE,
  grid = expand.grid(decay = 0)
)

evaluate(mlrf)

Multi-class area under the curve: 0.8178

Confusion Matrix and Statistics

          Reference
Prediction   AD ADE4  SCD SCDE4
     AD    5047 2100    4    91
     ADE4  3344  941   60   202
     SCD    115  215 1998  3271
     SCDE4   94  144 1938  5136

Overall Statistics
                                         
               Accuracy : 0.5313         
                 95% CI : (0.525, 0.5375)
    No Information Rate : 0.3522         
    P-Value [Acc > NIR] : < 2.2e-16      
                                         
                  Kappa : 0.3593         
                                         
 Mcnemar's Test P-Value : < 2.2e-16      

Statistics by Class:

                     Class: AD Class: ADE4 Class: SCD Class: SCDE4
Sensitivity             0.5869      0.2768    0.49950       0.5903
Specificity             0.8637      0.8307    0.82604       0.8640
Pos Pred Value          0.6969      0.2069    0.35685       0.7024
Neg Pred Value          0.7965      0.8780    0.89519       0.7950
Precision               0.6969      0.2069    0.35685       0.7024
Recall                  0.5869      0.2768    0.49950       0.5903
F1                      0.6372      0.2368    0.41629       0.6415
Prevalence              0.3482      0.1377    0.16194       0.3522
Detection Rate          0.2043      0.0381    0.08089       0.2079
Detection Prevalence    0.2932      0.1841    0.22668       0.2960
Balanced Accuracy       0.7253      0.5537    0.66277       0.7272

Decision Tree

tree <- multifit(
  X = X,
  y = y,
  model = "rpart2",
  ctrl = ctrl,
  grid = expand.grid(maxdepth=3)
)

evaluate(tree)

Multi-class area under the curve: 0.8176

Confusion Matrix and Statistics

          Reference
Prediction   AD ADE4  SCD SCDE4
     AD    5578 1904   56    95
     ADE4  2692 1241   58   155
     SCD    182  153 2704  5420
     SCDE4  148  102 1182  3030

Overall Statistics
                                         
               Accuracy : 0.5082         
                 95% CI : (0.502, 0.5145)
    No Information Rate : 0.3522         
    P-Value [Acc > NIR] : < 2.2e-16      
                                         
                  Kappa : 0.3445         
                                         
 Mcnemar's Test P-Value : < 2.2e-16      

Statistics by Class:

                     Class: AD Class: ADE4 Class: SCD Class: SCDE4
Sensitivity             0.6486     0.36500     0.6760       0.3483
Specificity             0.8724     0.86362     0.7220       0.9105
Pos Pred Value          0.7308     0.29932     0.3197       0.6791
Neg Pred Value          0.8229     0.89496     0.9202       0.7198
Precision               0.7308     0.29932     0.3197       0.6791
Recall                  0.6486     0.36500     0.6760       0.3483
F1                      0.6872     0.32892     0.4341       0.4604
Prevalence              0.3482     0.13765     0.1619       0.3522
Detection Rate          0.2258     0.05024     0.1095       0.1227
Detection Prevalence    0.3090     0.16785     0.3425       0.1806
Balanced Accuracy       0.7605     0.61431     0.6990       0.6294

XGBoost Forest

xgb <- multifit(
  X = X,
  y = y,
  model = "xgbTree",
  ctrl = ctrl,
  grid = xgb.grid
)

evaluate(xgb)

Multi-class area under the curve: 0.8355

Confusion Matrix and Statistics

          Reference
Prediction   AD ADE4  SCD SCDE4
     AD    6432 2000   19     3
     ADE4  2079 1263   35   184
     SCD     24    6 1559  3086
     SCDE4   65  131 2387  5427

Overall Statistics
                                          
               Accuracy : 0.5944          
                 95% CI : (0.5882, 0.6005)
    No Information Rate : 0.3522          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.4336          
                                          
 Mcnemar's Test P-Value : < 2.2e-16       

Statistics by Class:

                     Class: AD Class: ADE4 Class: SCD Class: SCDE4
Sensitivity             0.7479     0.37147    0.38975       0.6238
Specificity             0.8744     0.89211    0.84947       0.8386
Pos Pred Value          0.7608     0.35468    0.33348       0.6775
Neg Pred Value          0.8666     0.89891    0.87810       0.8039
Precision               0.7608     0.35468    0.33348       0.6775
Recall                  0.7479     0.37147    0.38975       0.6238
F1                      0.7543     0.36288    0.35942       0.6496
Prevalence              0.3482     0.13765    0.16194       0.3522
Detection Rate          0.2604     0.05113    0.06312       0.2197
Detection Prevalence    0.3423     0.14417    0.18927       0.3243
Balanced Accuracy       0.8112     0.63179    0.61961       0.7312

Model Comparison

Observations:

1. Adding serum metabolite information (either the full 230-metabolite matrix or its 6-factor projection) seems to increase the discriminatory power of the models.

2. Fitting 6 ML-estimated factors obtained by the FMradio package (cummulatively explaining 30% of variace) yields increased classification performance, serving as a valuable dimension reduction technique for high-dimensional data.

3. Looking at the confusion matrix and individual ROC curves, all models were able to discriminate better among certain classes (AD+E4/SCD+E4, AD+E4/SCD, AD-E4/SCD+E4 and AD-E4/SCD-E4) compared to others (AD+E4/AD-E4 and SCD+E4/SCD-E4).

aucs <- data.frame(AUC = c("Clinical features only" = benchmark$roc$auc, "Clinical features + 230 metabolites" = mlr$roc$auc, "Clinical features + 6 latent factors" = mlrf$roc$auc, "Decision Tree" = tree$roc$auc, "XGBoost" = xgb$roc$auc))

knitr::kable(aucs)

Network statistics

# # Plot the densities of centrality degree scores
# plot(density(DegreesAD1[, 2]),
#   col = "blue", xlim = c(-1, 8), xlab = "Degree", main = ""
# )
# lines(density(DegreesAD2[, 2]),
#   col = "red"
# )
# legenda <- c("AD class 1", "AD class 2")
# legend(5, 0.5,
#   legend = legenda, 
#   lwd = rep(1, 2), lty = rep(1, 2), col = c("blue", "red"), cex = 0.7
# )

Wilcoxon Signed Rank test between the network statistics of ApoE4 carriers and non-carriers. The test shows distinct centrality degrees and average number of negative and positive edges.

# Perform a Wilcoxon signed rank test

w <- furrr::future_map2(NetStatsE4yes[, 1:7], NetStatsE4no[, 1:7], ~ {
  wilcox.test(.x, .y, paired = TRUE, alternative = "greater")
})
p.values <- list(1:length(w))
for(i in 1:length(w)){
  p.values[i] <- w[[i]][['p.value']]
}
names(p.values) <- names(w)
p.values[p.values<0.05]

$PCE4YES.degree
[1] 0.001016491

$PCE4YES.nNeg
[1] 0.007252307

$PCE4YES.nPos
[1] 0.02813907

Distinct metabolite communities among ApoE4 carriers/non carriers

# Get the communities per class
set.seed(141516)
CommC1 <- Communities(PCE4NO,
  Vcolor = Colors,
  Vsize = 10, Vcex = 0.5, main = "ApoE4 non-carriers"
)

CommC2 <- Communities(PCE4YES,
  Vcolor = Colors,
  Vsize = 10, Vcex = 0.5, main = "ApoE4 carriers"
)